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Abstract 

We present a high-order compact finite difference approach for a class of parabolic partial 
differential equations with time and space dependent coefficients as well as with mixed second- 
order derivative terms in n spatial dimensions. Problems of this type arise frequently in 
computational fluid dynamics and computational finance. We derive general conditions on 
the coefficients which allow us to obtain a high-order compact scheme which is fourth-order 
accurate in space and second-order accurate in time. Moreover, we perform a thorough von 
Neumann stability analysis of the Cauchy problem in two and three spatial dimensions for 
vanishing mixed derivative terms, and also give partial results for the general case. The 
results suggest unconditional stability of the scheme. As an application example we consider 
the pricing of European Power Put Options in the multidimensional Black-Scholes model 
for two and three underlying assets. Due to the low regularity of typical initial conditions 
we employ the smoothing operators of Kreiss et al. to ensure high-order convergence of the 
approximations of the smoothed problem to the true solution. 


1 Introduction 

In the last decades, starting from early efforts of Gupta et al. [nmn], high-order compact finite 
difference schemes were proposed for the numerical approximation of solutions to elliptic [B 
H] and parabolic [B [H] partial differential equations. These schemes are able to exploit the 
smoothness of solutions to such problems and allow to achieve high-order numerical convergence 
rates (typically strictly larger than two in the spatial discretisation parameter) while generally 
having good stability properties. Compared to finite element approaches the high-order compact 
schemes are parsimonious and memory-efficient to implement and hence prove to be a viable 
alternative if the complexity of the computational domain is not an issue. It would be possible to 
achieve higher-order approximations also by increasing the computational stencil but this leads 
to increased bandwidth of the discretisation matrices and complicates formulations of boundary 
conditions. Moreover, such approaches sometimes suffer from restrictive stability conditions and 
spurious numerical oscillations. These problems do not arise when using a compact stencil. 

Although applied successfully to many important applications, e.g. in computational fluid 
dynamics [laiiiiiisiiB] and computational finance [3111121 HU], an even wider breakthrough of the 
high-order compact methodology has been hampered by the algebraic complexity that is inherent 
to this approach. The derivation of high-order compact schemes is algebraically demanding, hence 
these schemes are often taylor-made for a specific application or a rather smaller class of problems 
(with some notable exceptions as, for example Tele’s paper [14] 1. The algebraic complexity is 
even higher in the numerical stability analysis of these schemes. Unlike for standard second- 
order schemes, the established stability notions imply formidable algebraic problems for high-order 
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compact schemes. As a result, there are relatively few stability results for high-order compact 
schemes in the literature. This is even more pronounced in higher spatial dimension, as most of 
the existing studies with analytical stability results for high-order compact schemes are limited to 
a one-dimensional setting. 

Most works focus on the isotropic case where the main part of the differential operator is given 
by the Laplacian. Another layer of complexity is added when the anisotropic case is considered 
and mixed second-order derivative terms are present in the operator. Few works on high-order 
compact schemes address this problem, and either study constant coefficient problems [7] or specific 
equations [5]. 

Consequently, our aim in the present paper is to establish a high-order compact methodology 
for a class of parabolic partial differential equations with time and space dependent coefficients and 
mixed second-order derivative terms in arbitrary spatial dimension. We derive general conditions 
on the coefficients which allow to obtain a high-order compact scheme which is fourth-order accu¬ 
rate in space and second-order accurate in time. Moreover, we perform a von Neumann stability 
analysis of the Cauchy problem in two and three spatial dimensions for vanishing mixed derivative 
terms, and also give partial results for the general case. As an application example we consider the 
pricing of European Power Put Basket options with two and three underlying assets in the mul¬ 
tidimensional Black-Scholes model. The partial differential equation features second-order mixed 
derivative terms and, as an additional difficulty, is supplemented by an initial condition with low 
regularity. We use the smoothing operators of Kreiss et al. [T3] to restore high-order convergence. 

The rest of this paper is organised as follows. In the next section, we state the general parabolic 
partial differential equation in n spatial dimensions and give the central difference approximation 
for the associated elliptic problem. We then derive auxiliary relations for the higher-order deriva¬ 
tives appearing in the truncation error of the central difference approximation in Section [S] In 
Section m we give conditions on the coefficients of the partial differential equation under which a 
high-order compact scheme is obtainable. Semi-discrete high-order compact schemes in n = 2 and 
n = 3 space dimensions are derived in Section 0 Section |5] discusses the time discretisation. A 
thorough von Neumann stability analysis of the Cauchy problem in n = 2 and n = 3 space dimen¬ 
sions is performed in Section [T] In Section [8] we apply the schemes to option pricing problems for 
European Basket Power Put options and report results of our numerical experiments in Section 121 
Section [TO] concludes. 


2 Parabolic problem and its central difference approxima¬ 
tion 


We consider the following parabolic partial differential equation with mixed derivative terms in n 
spatial dimensions for u = u{xi, ..., a;„, r), 

n 

+ X! 

ij'=l 
i<j 

with initial condition uq = u(xi,... Xn,0) and suitable boundary conditions, with space- and 
time-dependent coefficients = ai(xi,... Xn,T) < 0, bij = bij(xi,... Xn,T), ct = Ci(xi,... Xn, r) 
and g = g(xi,... Xn, t). The spatial domain LI C M" is of n-dimensional rectangular shape with 
LI = III X ... X Ll„ and Xi G Lli = a:mL] with < a:mL for i e {1,..., n}. The temporal 

domain is given by LIt = ]0,Ti„ax] with Tmax > 0. The functions a(-,T), 6 (',t), c(-,t) and gi-,T) 
are assumed to be in C‘^{Ll) for any r G Llr, G C^{Ll) and u is assumed to be differentiable 

with respect to r. Introducing / := —Ur g we can rewrite © as 




XidXn 


E< 


du 
'' dxi 


=g in LI X Llr, 


( 1 ) 


+E- 


d'^u 

"dxj 


n 


n 


l—L = l 

i<j 


d'^t 


dxidxj 


ft 

1—1 


( 2 ) 
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We start by defining a grid on Lt, 

e n I x['^^ = xl^l^ + ikAxk,0 <ik < Nk, k = l,2 ,... ,n}, (3) 

where Axk = (a^max — xl^}^)/Nk > 0 are the step sizes in the fc-th direction with iVfc G N for 
k = l,2,...,n. We use for the interior of G^"^. On this grid we denote by the 

discrete approximation of the continuous solution u at the point G G^"^ and 

time T € Or- Using the central difference operator and the standard second-order central 
difference operator in ccfe-direction we get 


d'^u 

du 

dxk 

d‘^u 


=Dlu - 




{Axk)"^ d^u , 4 . 

—Di.u „ ^^3 -I- O {{Axk) ) 


^ ^ =DlDlu - 

dxkdxp ^ 


6 

(Axfc) 


k 

2 


(Aa 


d'^i 


6 


dxldxp 


6 dxkdxl 


0{{Axkf) 


+ O {{Axkf{Axpf) + O {{Axpf) + O 


/ {Axkf 

\ Axp 


(4) 


for k,p G {1, 2,..., n} and k ^ p, evaluated at the grid points 3 :-^^ ..., G G^"\ Using 
the approximations (U) in @ gives 


n n n n 

/ = ^ a^D^u + ^ ^ ^ 

i<j 


ai{AxiY 
12 ^ 


i ,a=i 

i<j 


(Axi)'^ d‘^u 
6 dxfdxj 


(Axj)'^ d'^u 
6 dxidxj 


-E 


Ci(Aa:i)^ d^u 

6 dx^ 


+ s, 


(5) 


where e G O (/i^) if Axi G O (h) ior i = 1,2,... , n for a step size h > 0. If the consistency error 
is in O (h'^) , we call the scheme high-order. In order to achieve a high-order scheme we need to 
find second-order approximations of the derivatives and . for *, j G {1,..., n} with 

i ^ j- We call the scheme high-order compact, if we can achieve this using only points from a 
compact computational stencil for x = (a;-^\ a;-^\ ..., G G^"^ We have 

U{x) = { ..., a:,Efc„) e G(”) I G {-1, 0,1} for m = 1, 2,..., n} (6) 

for x = ..., x^^^) as the compact computational stencil and define ~ u(x^^\ x^^^ ..., a;^"^). 


3 Auxiliary relations for higher derivatives 


In this section we calculate auxiliary relations for the higher derivatives appearing in ([S])- These 
relations for the higher derivatives can be calculated by differentiating @. In doing so no addi¬ 
tional error is introduced. Differentiating equation m with respect to Xk and then solving for 
leads to 


d^u ^ Gi d^u ^ 1 dai d'^u 1 dau d'^u ^ hij d^u 
dxl Gk dxfdxk A;' Ufe dxk dxj Gk dxk dxl Gk dxidxjdxk 

i^k i^k i<^j 


^ 1 dbij d'^u 

Gk dxk dxidxj 
1,3 = 1 ■> 

i<3 


E 


Ci d'^u 
Gk dxidxk 


E 


1 dci du 
Gk dxk dxi 


1 9f 

— =-^ 
Gk OXk 


k 


(7) 
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for k = 1,... ,n. The relation for Ak can be approximated with consistency order two on the 
compact stencil (|6]) using the central difference operator, as all derivatives of u in the above 
equation are only differentiated up to twice in each direction. 




Differentiating m twice with respect to Xk, and solving the resulting equation for 
obtain 


we 


d^u 

d4 


= -E 


2=1 

i^k 


d^u 


Qk dxfdxl 
1 d'^u 


+ 


2 dtti d^u 
ttk dxk dxfdxk 


Uk dxl dxl 


k-l 


-E 




d‘^1 


1 d'^u 

ak dxl 

2 dbii 


2 dak d^u 
ak dxk dxl 


d^i 


1 d'^bij d'^u 


ij=l 

i<j 

k-l 


ak dxidxjdxl ak dxk dxidxjdxk ak dxl dxidxj 


bik d^u 

^ ak dxidxl ^ 


i=l 


i=l 


2 dbik d^u 1 d'^bik d'^u 
ak dxk dxidxl ak dxl dxidxk 


- E 




d'^i 


t=fc+i 


ak dxjdxl 


-E 


d^i 


- E 

^ j—k+l 
2 dcj 


2 dbkj d^u 


+ 


1 d'^bkj 9^ 


ak dxk dxjdxl ak dxl dxjdxk 


dh 


. ttfc dxidxl 
' ^ flfc dxidxl 


^ 1 d'^Ci du 
ak dxk dxidxk ak dxl 

bkj d‘^u 


ak dxl 


- E 

i=fe+i 


Ofc dxjdxl 


( 8 ) 


We can approximate Bk with second order consistency on the compact stencil (jS]) , when using the 
central difference operator and the auxiliary relations for Ak in © for k = 1,... ,n. Differentiating 
equation m once with respect to Xk and once with respect to Xp leads to 


ak 


d'^u 

dxldxp 


= -E 


d^u 

'^dxkdxl 

d^u 


dai d^i 


dai d^u 


d^a, d^i 


i=l 
i^k,p 

dap d^u 


dxldxkdxp dxkdxldxp dxpdxldxk dxkdxp dxj 


dap d^u 
dxk dxl 


9^ Op 9 ^m dak d^i 


dxp dxldxk dxkdxp dxl 

d^u db. 


dxk dxldxp 


dak d^u d'^ak d'^u 
dxkdxp dxl 


-E 


+ 


dh 


+ 


dxp dxl 


dbi 


d^i 


d%,, d^i 


i,j=i 

i<j 


dxidxjdxkdxp dxk dxidxjdxp dxp dxidxjdxk dxkdxp dxidx 


-E 


d^i 


dci 9^z 


9c, 9^1 


9^C,; du 


' dxidxkdxp dxk dxidxp dxpdxidxk dxkdxpdxi 


+ 


9V 

dxkdXf 


=:C> 


kp') 


where Ckp can be approximated on the compact stencil ([HI) using Ak and Ap, as defined in equation 
and the central difference operator for k,p = 1,..., n with k ^ p. This can be written as 


d\ 


dxldxp 


Ckp 

ak 


d'^u 


ak dxkdxl' 


(9) 


4 Conditions for obtaining a high-order compact scheme 

In this section we derive conditions on the coefficients of the partial differential equation © 
under which it is possible to obtain a high-order compact scheme, i.e. only using points of the 
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n-dimensional compact stencil m for discretisation and receiving a fourth-order scheme with 
Axi S O (/i) for j = 1,..., n for a given step size h > 0. Using equations d?]) and dH) and then ([9]) 
in m leads to 


n n n "/"AUR 

/ = ^ UiDju -f ^ -I ^ CiD'^u - ^ ^e 


i,j = i 

i<j 


12 


-E 


biji^Axi^ Ci^ 


jj'=l 

i<j 


12 a,; 


-E 


bi, d‘^1 


12 dxidx^ 

i<j 


[Ax.r - 


2 aj{Axi) 


21 


-E 


2=1 


c^{Ax^)^A, 

6 ’ 


( 10 ) 


where e G 0(h‘^), if Axi € O (h) for i = 1,... , n. The leading error terms are given by 
b., a u [Axj)^ — Jqj. j g { 1 ^... ^ n] with i ^ j. If the conditions 


12 dxidx^ L 


bij =0 or 


(Axi )2 


^3 


( 11 ) 


are fulfilled for all i, j G { 1 ,..., n} with i ^ j these second order terms vanish and the resulting 
error term is of fourth order. Hence, for any partial differential equation m which satisfies m 
we obtain a high-order compact scheme. In the case bij = 0 for all i,j G 1,... ,n, it is possible 
to choose Axi > 0 freely for each spatial direction, whereas in other possible cases there are 
interdependencies for at least some of the step sizes. For each pair {i,j) with bij ^ 0 the condition 

^ has to hold for all x = ..., G G^'^\ This means aj/ai has to be 

constant as (Aa;j)^/(Aa:i)^ is constant, see dS])- 


5 Semi-discrete high-order compact schemes 

In this section we present the semi-discrete high-order compact schemes in spatial dimensions 
n = 2,3. We consider the case where the cross derivatives do not vanish, hence we assume, for 
simplicity, a^ = a in combination with Axi = h > 0 for i = 1 ,... n to satisfy condition dm. Our 
aim in this section is to derive a semi-discrete scheme of the form 

^ K^{x,T)Ut^^...^t^{T)] =g{x,T) (12) 

xeG('*) 

° (n) . 

at each point x G G with Ax* = h > 0 for i = and time r, where the fnnction 

° (n) 

g : G X ^ R depends on the function g given in dH). 

5.1 Semi-discrete two-dimensional scheme 

In this section we derive the high-order compact discretisation of m in spatial dimension n = 2 . 

Considering the grid point G with Axi = Ax 2 = h > 0 and time r G we are 

able to obtain the coefficients Ki^m of Ui^m (t) for 1 G {ii — 1, ii, ii + 1} and m G {*2 — 1, * 2 , ^2 + 1} 
on the compact stencil by employing the central difference operator in (HOI) . To streamline notation 
we denote by [-Jit the first derivative with respect to Xk and by the second derivative, once in 
Xk- and once in Xp-direction with k,p G {1, 2}. Note that in the following the functions a, 61 ^ 2 , ci, 

C 2 and g are evaluated at G and r G U,-. We omit these arguments for the sake of 
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readability. The coefficients are given by: 


Eii .io — 
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Analogously, we obtain the coefficients Mi,m of drill,m (t) for / € {ii — + 1} and 

{*2 — + 1 } at each point G and time r G Or, 

M -M M 1 _ h[a]2 _ bi2h[a]i , cah 

1 ’2 12 ^ 24a2 24a 12 a ^’2 3 

^ ( 2 ) 

where Axi = Ax 2 = h > 0. Additionally, for a; G G , t G Or, 

(h'^a^ci-2h'^a'^[a]i - bi 2 h'^[a] 2 a) [g]i h‘^[g]ii bi 2 h‘^[g]i 2 
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holds, where Axi = Ax 2 = h > 0 was used. We have , r) = Kn^ ,n 2 and {xni , Xn}, t) 

Mni,n 2 in (HD) with ni e {ii - l,iiAi + 1} and n 2 G {12 - + 1} for x = {xf^,xf^) G 

and T G Ot- Kx and are zero otherwise and the approximation only uses points of the compact 
stencil. 


5.2 Semi-discrete three-dimensional scheme 

In this section we derive the high-order compact discretisation of o in spatial dimension n = 3. 
Considering the conditions in dill) we observe that in the three-dimensional case we have three 
different possibilities to satisfy the conditions and thus obtain a high-order compact scheme. We 
focus on the case a = oi = 02 = 03 and set h = Axi = Ax 2 = Ax^. Considering an interior grid 

point G and time t G Ox we are able to produce the coefficients Kk,i,m of 

Uk,i,m (t) for /c G {zi - 1, ii, h + 1}, ^ G {z 2 - 1, 12 , *2 + 1} and m G {13 - 1, 13,13 -f 1} by employing 
the central difference operator in (fTUl) . Again, to streamline the notation we denote by [-Jfe and 
[•]fcp the first and second derivative of the coefficients with respect to Xfc, and with respect to Xk 
and Xp, respectively. Note again that in the following a, & 12 , 613 , & 23 j ci, C 2 , C 3 and g are evaluated 

at (x|^\ x-^\ x-^^) G and t G Ox, where Axi = > 0 for z = 1, 2,3. We omit these arguments 

for the sake of readability. Due to the length of the coefficient expressions Kk,i,m, they are given 
in the appendix. 

In a similar way we define Mkj^m as the coefficient of dxUk,i,m (x) for 
fc G {zi - l,zi,zi -I- 1}, Z G {z 2 - 1 ,Z 2 ,Z 2 + 1} and m G {is - 1 ,Z 3 ,Z 3 + 1} by 


.^n±i,z 2 -i,i 3 ~-^nTi>» 2 +i,i 3 “ ^ 48 a’ ~ 2’ 

^ ^ ^13 ^ ^ ^23 

fo^n±l,Z 2 ,j 3 -l “fo^ilTl.» 2 ,i 3 + l “ ^ 4 ^’ * 2 ± 1 , 13-1 “ fo^il,j 2 Tl>i 3 ±l “ ^ 4 ^’ 

4 h\„^. 

fo^ii±l,i 2 .i 3 Y2 ^ 

-^ii,i 2 ±l,i 3 “42 ^ 

-^il,j2,j3±l T 


24a2 

24a2 

± -=F —^ 

24a 12a 

^^12 [a] 1 
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^ hc2 h[a]2 

24a2 
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^^23 [a] 2 

hbi3[a\i 

^ hc3 ^ h[a\3 


fo^ii±l.i2-l,j3-l —fo^ii±l,i2±l.Z3-l — fo^ii±l.i2-l,i3±l ” ±l,i2±l,i3±l ” 0. 


For the right hand side of (IT^ we have for x = (x-J\ x^^\x-^^) G t G Ox, 

_ icih?‘a-2h‘^[a]ia-bi2h?[a]2-bish?‘[a]s)[g]i bishT^Mis 

{c2h^a - 2h‘^[a\2a - &i2/z^[a]i - b23h'^[a]3) [gh &23^^[5]23 
12^2 + 

[csK^a - 2h'^[a]sa - hish'^[a]i - 623 / 1 ^H 2 ) [ 5)3 h‘^[g]ii 
^ m 2 + 42 

hl 2 h?'[g]l 2 fe^[g]33 fe^[g ]22 

12a 12 12 

We define Kx{xn} ,Xn} ,Xn} ,t) = Km,n 2 ,n 3 and Mx{xn} ,Xn} ,Xn} ,t) = for each point 

X = {x[^\ x['^\ x[^'^) G G and r G Ox, where rij G {ij — + 1} with j = 1,2,3. and 

Mx are zero otherwise. Hence, the approximation only uses points of the compact stencil ®. 
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6 Fully discrete scheme 

The semi-discrete scheme presented in the previous sections can be extended to a fully discrete 
scheme for the parabolic problem m by additionally discretising in time. Any time integrator 
can be implemented to solve the problem as in [^D]. Here we consider a Crank-Nicolson type 
time-discretisation with constant time step At to obtain a fully discrete scheme. Let 

'' At '' At 

A^{x,Tk+l) = {x,Tk) -L —Kx {x,Tk+l) , Bx{x,Tk) = Mx {x,Tk) - {x,Tk) , 

where Mx{x,Tk) = {Mx {x,Tk) -I Mx {x,Tk+i)) /2. Kx{x,t) and Mx{x,t) are defined through 
a semi-discrete finite difference scheme with fourth-order consistency using only points of the 
compact stencil ©• Then, a fully discrete high-order compact finite difference scheme for © with 
n G N on the time grid Tk = kAr for fc = 0,..., Nr and Axi = h for all i is given at each point 

X = e by 

Ax{x,Tk+i)Uf{+{i^= Y Bx{x,Tk)Ull^^^^^i^ + ^g{x,Tk,Tk+i), ( 13 ) 

xg(J(x) xGU(x) 

where g {x,Tk,Tk+i) = g{x,Tk) -\-g{x,Tk+i) and x = ... , G U (x). This scheme is 

second-order consistent in time and fourth-order consistent in space. We have fourth-order con¬ 
sistency in terms of h for At G 0{h?) while only using the compact stencil. Note that up to this 
point only the spatial interior is discussed. The applied boundary conditions may still have an 
effect the above numerical scheme. 


7 Stability analysis for the Cauchy problem in dimensions 
n = 2,3 


In this section we consider the stability analysis of the high-order compact scheme for the Cauchy 
problem associated with © in the case n = 2,3. The coefficients of the semi-discrete scheme are 
given in Section 15.11 for two spatial dimensions and in Section 15.21 when three spatial dimensions 
occur. Those coefficients are non-constant, as the coefficients of the parabolic partial differential 
equation © are non-constant. 

We consider a von Neumann stability analysis. Other approaches which take into account 
boundary conditions like normal mode analysis m are beyond the scope of the present paper. 
For both n = 2 and n = 3, we give a proof of stability in the case of vanishing cross derivative terms 
and frozen coefficients in time and space, which means that all possible values for the coefficients 
are considered, but as constants, hence the derivatives of the coefficients of the partial differential 
equation appearing in the discrete schemes are set to zero. This approach has been used as well in 
[IIIIII] and gives a necessary stability condition, whereas slightly stronger conditions are sufficient 
to ensure overall stability m This approach is extensively used in the literature and yields good 
criteria on the robustness of the scheme. In (IT^ we use 

n 

with Sn=Y 

m—1 

for jm G {im ~ 1) *m) *m +1}) where I is the imaginary unit, g^ is the amplitude at time level k and 
Zm = 2'Khf\m for the wavelength Am G [0, 27r[ for m = 1,..., n. Then the fully discrete scheme 
satisfies the necessary von Neumann stability condition for all Zi, Z 2 , when the amplification factor 
G = g^^^ jg^ satisfies 


|Gp - 1 < 0. 


( 14 ) 
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7.1 Stability analysis for the two-dimensional case 

In this section we perform the von Neumann stability analysis for the two-dimensional high-order 
compact scheme of Section [53] The analysis of the case with vanishing cross-derivative and frozen 
coefficients are carried out in detail. In the case of non-vanishing mixed derivatives partial results 
are given for frozen coefficients. 

Theorem 1. For a = ai = 02 < 0, 61^2 = 0 and Axi = Aa ;2 = h > 0, the fully discrete high-order 
compact finite difference scheme given in (1131) with n = 2, with coefficients defined in Section \5.1[ 
satisfies (for frozen coefficients) the necessary stability condition (1141) . 

Proof. Let = cos(2i/2), ^2 = cos(z2/2), pi = sin{zi/2) and r ?2 = sin(z2/2). The stability 
condition (HI for the fully discrete scheme m using the coefficients defined in Section [53] yields 
[Gp — I = Ng/Dq (explicit expressions for Nq, Dq are given below). We discuss the numerator 
Ng and the denominator Dg separately in the following. 

The numerator can be written as Ng = Ska (jiih"^ -\- n 2 h'^) where the polynomials 

712 =8a^/i (Cl) '? 2 ) /2 (Cl) C 2 ) and 774 = /a (Ci) /4 (Ci) C 2 ) c? -f /a (C 2 ) A (C 2 ) Ci) 


are non-negative, since 


/2 {x, y) =2 - -f - y > 0, 

fi {x, y) =2x^y‘^ - - 1 < 0 , 


fi {x,y) =x^ +y^ + 1 > 0, 

/a (x) =x^ - 1 < 0, 

for x, 2 / G [—1,1]. Hence, we observe that Ng < 0 holds, as Ci)C 2 S [—1) !]• 

Now we consider the denominator Dg, which can be written as 

Dg = deh^ -\- ((i4,2^^ + 1 ^ 4 , 1 ^ + 1 ^ 4 , 0 ) T ( 1 ^ 2 , 2 ^^ T d2^ik)h’^ -|- dp, 

where 

do =16a^k^ (2C?Cl + C? + Cl - 4)" > 0, d 2 ,i = 16a3/i (^i, 6) /s (Ci) C 2 ) > 0, 
d2,2 =4a^ 9 (Cidici -I- C2d2C2)^ + 2/a (Ci) /e (Ci)C2) c| -I- 2/a (C 2 ) /e (C2,Ci) c| ) 
d4,o =4a^/i (Ci)C2)^ > 0, d4,i = -4an4 > 0, 

d4,2 = [/3(Cl)c| — 27yi772ClC2ClC2 + f3if.2)cl] >0, dp = (CldlCl + C2d2C2)^ > 0, 
because a < 0 and where 


h {x, y) =2x^y^ -f x^ -f - 4 < 0, /p (x, y) = 2x^y'^ - 5x^ - y^ -f 4 


with x,y G [—1,1]. We observe that fe{x,y) changes sign, as, for example /s(0,0) = 4 and 
/p (1, 0) = —1. Hence, we cannot determine the sign of d 2,2 directly. 

If Cl = C 2 = 0, we have d 2,2 = 0 and hence Dg > 0. Since d 2,2 is symmetric, we can say 
without loss of generality that ci ^ 0 in the following. Furthermore, as both ci and C 2 are frozen 
coefficients, we set m = C 2 /C 1 , which leads to 

d 2,2 =4a^Ci[9(Ci7?i +f 2 V 2 mf + 2/a(Ci)/p(Ci)C 2 ) + 2/a(C2)/6(C2,Ci)"i^] =: la^clgijn). 


The function g (m) can be rewritten as 

9 {m) =vlf7 (Cl, C 2 ) + 18CiC2did2TO + vlfr (C 2 , Ci) 

with /y (x, y) = 4x‘*y^ — 2x^ — y^ + 8 > —2x^ — y^ -|- 8 > 5. In the case yi = 0 we have y(m) = 
dl /7 (C 2 , Cl) > 0 and thus d 2,2 > 0 and Dg > 0 . In the case yi 7 ^ 0 we have y|/ 7 (Ci, C 2 ) > 0 , hence 
the function g (m) has a global minimum. This minimum is located at 


—9ClC2d2 

dl/7(Cl)C2)’ 


which leads to g (rh) 


2dl/5 (Cl)C2) /s 

fr (Ci)C 2 ) 
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where /g = + C? + Cl “ 2CfCl^l “ - 8 < 0. Since /s (^ 1 , 6 ) < 0 we have g(m) > 

0 for all m G K, and thus we have Dc > 0 for all cases as a < 0 . 

We still need to show that Dq > 0 for all Ci,C 2 G H holds do > 0 for all (^ 1 ,^ 2 ) G 

[—1,1]^ \ { — 1,1}^ as a < 0 and k > 0. This leads to Dq > 0 in these cases. For the case 
(Ci,C 2 ) G { — 1, Ip it holds /i (Ci,C 2 ) = 3, which leads to ^ 4^0 = 36a^ > 0 and Dq > 0. Therefore, 
we have Dq > 0 for all (^ 1 ,^ 2 ) G [—1,1]^ and condition (ITTl) is satisfied. □ 

For 61^2 ^ 0 the situation becomes much more involved. Many additional terms appear in the 
expression for the amplification factor G and we face an additional degree of freedom through 
&i^ 2 - Since we have proven condition (HI holds for 6 i _2 = 0 it seems reasonable to assume it also 
holds at least for values of 61.2 close to zero. In von Neumann stability analysis it is often most 
difficult to guarantee that stability condition (O holds for extreme values of ryi, 772 , Ci and ^ 2 - We 
have the following partial result which holds in the case of frozen coefficients and non-vanishing 
coefficient of the mixed derivative, i.e. 61^2 ^ 0 . 

Lemma 2. For a = oi = 02 < 0, arbitrary 64^2 o,nd Aa:i = Aa ;2 = h > 0, the high-order compact 
scheme (US with the coefficients for the two-dimensional case defined in Section \5. 1\ satisfies (for 
frozen coefficients) the stability condition (1141) at the corner points Ci = il and C 2 = il- 

Proof. Using 771 = sin(zi/2) = — Ci = 0 for Ci = ±1 and 772 = sin(z 2 / 2 ) = \/l — Cl = 0 

for C 2 = ±1, straight-forward computation shows that on each corner point |Gp — 1 = 0. Hence, 
condition (ITTl) holds. □ 

It is worth mentioning that in a comparable situation in (where a specific partial differential 
equation of type ([T]) is considered) an additional numerical evaluation of condition (IT4l) revealed 
it to hold also for non-vanishing mixed derivatives with (CiiCl) 7 ^ (1) !)■ However, the left hand 
side of m was very close to zero, and although the inequality was always satisfied, this left little 
room for analytical estimates. This leads to the conjecture that the stability condition in that 
case was satisfied also for general parameters, although it would be hard to prove analytically. 
Lemma[2]above suggests the present case is similar. We remark that in our numerical experiments 
we observe a stable behaviour throughout, also for general choice of parameters. 

7.2 Stability analysis for the three-dimensional case 

In this section we analyse the stability of the high-order compact scheme with coefficients given in 
Section r5.2l in three space dimensions. We first perform a thorough von Neumann stability analysis 
in the case of vanishing cross derivative terms for frozen coefficients. We observe no additional 
stability condition in this case. Then we give partial results in the case of non-vanishing cross¬ 
derivative terms for frozen coefficients. 

Theorem 3. For Oi = a < 0, bij = 0 and Axi = h > 0 for i,j G {1, 2,3}, i yf j, the fully discrete 
high-order compact scheme given in (na) with n = 3, with coefficients given in Section [KM satisfies 
(for frozen coefficients) the necessary stability condition (|TT)) . 

Proof. Let = cos{zi/2) and iji = sin(zi/2) for i = 1,2,3. The stability condition (|T4)) yields 
|Gp — 1 = Ng/Dg (explicit expressions for TVg, Dg are given below). 

For the numerator we have Ng = —Safe ( 774 / 1 "^ -I- 772 / 1 ^) < 0, since a < 0 and the polynomials 

772 =4a^/l (Cl) 6 ,^ 3 ) [/2 (Cl)C 2 ) + /2 (C3)Cl) + /2 (^2,^3)] < 0 , 

^4 = [/3 (6)6) + h (6)6)] c\ + [/a (6)6) + h (6)6)] c| + [fs (6)6) + 6 (6)6)] c| 

- Vs ( 6 »?iCi + 6 ?? 2 C 2 )^ - vl ( 6 ??iCi + fsVsCsf - Vi ( 6 ?? 2 C 2 + fsVsCsf < 0, 
are non-negative since 

/i {x, y) =x^ + + 2 :^ > 0, /2 [x, y) = 2x^7/^ - - 1 < 0, 

fs (x, y) =x‘^y‘^ (1 - x"^) -\- (x^ - l) < (l - x^) -k (x^ - l) = 0, 
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for x,y,z e [-1,1]. 

The denominator Dq can be written as 

Dq = df)h^ + {di^2k^ + di^\k + c?4,o) T (^ 2 , 2 ^^ + d 2 ,ifc) h"^ + doi 

where 

do =16a'‘fc^ [toi(^i, 52 ) + wi(^ 3,5 i) + mi(^2,6)]^ > 0, ^2,1 = 4an2 > 0, 

d2.2 =4a^ [me (Ci,di>6)c? + 2m7 (6) Ci6did2CiC2 + me {^2,V2,^i)cl 
+ me (6,di.6)ci +2m7 (6) CiCsdidaCiCa + me (C3,d3,6)c3 
+ me (6.d2,6)ci + 2m7 (6) 66d2d3C2C3 + me {^3,V3,^2)cl 
Tms (?7i, 6 , 6 ) Cl + ms (7?2, 6 , 6 ) C 2 + ms (7?3, , $ 2 ) C 3 ] 

d 4 ,o = 4 a^m 2 ( 6 , 6 ) 6 )^ > 0 , ^4.1 = 4 on 4 > 0 , de = KidiCi + 6 d 2 C 2 + ^ 311303 ]^ > 0 , 
d 4,2 = [? 7 iCi + ? 72 C 2 + 77303 + 2 ^i? 7 i^ 2772 CiC 2 + 2 ^i? 7 i^377301C3 + 2527726^30203] ^ > 0 , 
since a < 0 and 


mi {x, y) =2x^y‘^ — x^ — 1 < x^ — 1 < 0, 111,2 {x, y, z) = x^ + y'^ P > 0, 
7773 (a;, y) =x^y'^ (l - x^) + 7 /^ {x^ - l) < 7 /^ (l - x^) + 7 /^ {x^ - l) = 0, 
m4 (x, y) =(1 - x^)[x^{y^ - 1) + y^{x^ - 1)] < 0, 
ms {x, y,z) =- Sx'^y'^z'^ + Ix^y'^z^ + 4a;^ > -8x^7/^z^ + Ix'^y'^z'^ + Ix"^ 

= - 4a;^7/^z^ + 4x^ > -4x^ + 4x^ = 0, 

me (xi, X2, y) =4x3X1/ + (-8x3X1 + 2x3)/ +xl+ “^xlxl e [0,3], 
m7 (x) =2x^(x^ ~ (1 ~ x^)) + 7 > 0, 


for x,y,z G [—1,1]- We still need to show d2,2 > 0. Since we cannot determine the sign of d2,2 
directly, we consider three different cases. 

Having 5| = 51 = 1 leads to 

d2,2 =4a^ [2(-2.55 i77i + 8773)01 + (-8773 + 8773) C3] > 0 

as 5i < 1 and 773 < 1. 

Secondly, we consider ci = C2 = 03 = 0. This leads directly to d2,2 = 0. 

From now on we have (01,03,03) ^ (0,0,0). Since d 2,2 is symmetric with respect to 01,03,03, 
we assume without loss of generality oi yf 0. Additionally, we have (52)^1) 7^ Setting 

P 2 := C2/01 and ps := C3/01 gives 

d2,2 =40^03 [me (5i,77i,52) + 2m7 {^ 3 ) ^i^2mV2P2 + me {^2,V2 ,^i)pI 

+ me (5i, 771,53) + 2m7 (52) 5i677i V3P3 + me (53,773, 5i) Ps 

+ me ( 52 ,772,53) pI + 2 m 7 ( 5 i) 5253 772773P2P3 + me ( 53 , V 3 ,b) pj 
+ms (771,52,53) + ms (772,5i, 53) P2 + ms (773, 5i, 52) P3] 

=:4a^Oi [A:iiP 2 + ^22^3 + ki2P2P3 + kiP2 + ^2^3 + fco] =: 4a^03g {p2,P3) ■ 

To calculate the extremum of g {p 2 ,P 3 ), 


Vg(p2,6) 


/ 2kiiP2 + ki2P3 + ki\ 
\k 12 P 2 + 2.k22P3 +k2) 


is necessary, which leads to 


2kik22 — k2ki2 

kh - 4fc?ifci2 


2^261 — kiki2 

kh - 4fci/i2 ’ 



P3 = 


where ^33 — 4^33^33 = 771 <72 <73 
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with 


9l 92 = —2 Cl^C2^ — 2 + ^2^ + ^3^ + 3 G [0,4], 

93 =8a%"6" +4a'e2'6" +46^6^ - 46^6" 

- 4 ^ 4 ^ 3 ^ - 22 - 6 - 6 6 ^ 6 " + 8 

+ 8 + 20 ^ 2 %" - 2 6 " - 3 6 ' - 3 Cs" - 6 G [- 9 , 0 ], 


It holds 919293 7^ 0 for (CliCl) (IjI)- Since this is the unique root of Vg, as kii,k22 > 0 , we 
have a minimum at (p2,Pz) = iP2,P3)- We obtain g {p2,P3) = 9495/961 where 



6 g [ 0 , 9 ], 

4c^6" - 22e,e2e3 


with 96 7^ 0 for ^ ( 1 > !)• Hence, in all three cases we conclude c?2,2 > 0, and Dq > 0 holds. 

We still need to show that Dq > 0 for all 6i?2,C3 G [“ 1 j !]• H holds do > 0 for all (6i?2,^3) G 
[— 1 , 1 ]^ \ {— 1 , 1 }^ as a < 0 and k > 0 . This leads to Dq > 0 in these cases. For the case 
(6,61^3) G {“Ij 1 }^ we have m2 (61^2,^3) = 3 , which leads to d4_o = 36 a^ > 0 and Dq > 0 . 
Therefore, Dq > 0 holds for all (6i?2,C3) ^ [~ 1 j 1 ]^ and condition ([TT)l is satished. □ 

For the more general case with non-vanishing cross-derivatives we have the following result. 
The comments made in the previous section also apply here. 

Lemma 4. For Oi = a < 0, Axi = h > 0 for i = 1,2,3 and arbitrary bi^ 2 , ^ 1,3 62 , 3 , 

the high-order eompact scheme (ESI) with the coefficients for the three-dimensional case defined 
in Section \5.S\ satisfies (for frozen coefficients) the stability condition (fTH) at the corner points 
6 = il, 6 = 4=1 and 6 = 4=1- 

Proof. Using sin (2:1/2) = 1/I — 6 = 0 for 6 = 4 = 1 , sin(z2/2) = ^1 — = 0 for 6 = 4=1 
and sin(z3/2) = y^l — = 0 for 6 = 4 = 1 , straight-forward computation yields just as in the 
two-dimensional spatial setting to \G\^ — 1 = 0 for all corner points. Hence, condition (|T 4 )) is 
satisfied. □ 


8 Application to Black-Scholes Basket options 

To illustrate the practicality of the proposed scheme we now consider the n-dimensional Black- 
Scholes option pricing PDF (see, e.g. [ 23 ]). In the option pricing problem mixed derivatives appear 
naturally from correlation of the underlying assets. After transformations, the conditions (HID are 
satisfied, and we give the coefficients of the resulting scheme. Then we discuss the boundary 
conditions as well as the time discretisation. 

8.1 Transformation of the n-dimensional Black-Scholes equation 

In the multidimensional Black Scholes model the asset prices follow a geometric Brownian motion, 

dS^{t) = {pi - Si)S^{t)dt -\- a^S^{t)dW^{t), ( 15 ) 

where Si is the Uth underlying asset which has an expected return of pi, a continuous dividend 
of 6 , and the volatility CTi for i = 1 ,..., n and n € N. The Wiener processes are correlated with 
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{dWi, dWj) =: pijdt ior i,j = with i ^ j. Application of Ito’s lemma and standard 

arbitrage arguments show that any option price V{S, a, t) solves the n-dimensional Black-Scholes 
partial differential equation, 


dt 


2 = 1 


d'^V 

dS? 


E 

*,1=1 

i<j 


d'^V 


E 

2=1 


ViSt 


dS, 


- rV =0, 


where r]i = r — Si. The transformations 


Xi =j\n{Si/K)/ai, T = T — t and u = e'''^V/K, 


(16) 


(17) 


for i = 1 ,..., n, where 7 is a constant scaling parameter to assure that the resulting computational 
domain does not get too large, leads to 

n 

- 7^ E 

j.j=i 
i<j 

where i;i = ai/2 — iji/ai. Comparing this equation with ([T]), we identify 

Oi = — b,j = -TPij, Ci = 7 Ci, = 0, (19) 

for i, j = 1,..., n and i < j. We find that the transformed partial differential equation m with 
these coefficients satisfies the conditions given by m, if Axi = h for a step size h > 0 is used. 
Hence, we are able to obtain a high-order compact scheme in any spatial dimension n G N. 

We consider a European Power-Put Basket option, thus the final condition for (ITBl) is given by 




dxidxj 


n r. 

E UU 

2=i 


= 0 , 


(18) 


Ur — 


7 V- 


d^i 


2 

2=1 ^ 


V{Si,...,Sn,T) =m.a.x(K-'^uJiSi,o\ 
^ 2=1 ^ 


n _ 

where p is an integer and the asset weights satisfy ^ Wi = 1. Applying the transformations (I17|) 

i=l 


leads to the initial condition 


/ " cr X 

u(xi, ■. ■, Xn, 0) =Ar^~^ max( 1 — y^cuje ^ ,0) . (20) 

i=l ^ 

8.2 Semi-discrete two-dimensional Black-Scholes equation 

In this section we apply our general two-dimensional semi-discrete scheme, see Section 15.11 to the 
two-dimensional Black-Scholes model. To obtain the semi-discrete scheme m we have to apply 
(HU) with n = 2 to the coefficients in Section [571 which gives 

7^(5-2pf2) _ 7^12 I 7^1 7‘?2Pi2 7^ 

3h2 3 ’ 3h2 3h^ 3h 6 3h^' 

7 VL ^ 7^ ^ 7‘?iPi2 7^ 

3^2 3h 3h 6 3/i2 ’ 

, S291 _ 2^ 7£L _ 7‘?1P12 , 7*^2/012 _ 7^ , 7^Pl2 _ 7^Pl2 

12 12h 12h 6h 6h 12^,2 4/i2 Qh'^ ’ 

7^2 ‘^2‘^1 , 7^1 '1P12<^1 I 7‘^2P12 7^ 7^Pl2 7 ^Pi2 

12h ^ 12 12h 6h 6h 12h^ ^ 4/i^ 6h‘^ ’ 


-^21,22 ~ 
-^2l,22=tl 
^21 it 1,22 — 1 
-^2llhl,22 + l 
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where is the coefficient of Ui^m (t) for I € {ii — 1, ii,ii + 1} and m G {i 2 — 1, * 2 , *2 + !}• The 
coefficients of drUi^m (t) are given by 


~ q ’ 

^n±i.*2 ^2 127’ 


-Tfii + l,i2±l —^ 24 ’ 
,, _ 1 /l<?2 

*^’*2=^1 “12 ^ 127' 


Additionally, it holds g{x,T) — 0. This gives a semi-discrete scheme of the form (TT^ . where 
and Mx are time-independent. As in [5] we apply Crank-Nicolson type time discretisation and 
obtain the fully discrete scheme for the spatial interior. 


8.3 Semi-discrete three-dimensional Black-Scholes equation 

In this section we give the semi-discrete scheme (HID for the three-dimensional Black-Scholes Basket 
option. Using (TT^ with n = 3 in Section IST] and the appendix we obtain the coefficients Kk,i,m 
of Uk,i,m (t) for k G {ii — 1, ii, ii + 1}, Z G {12 — 1, *2, *2 + 1} and m G {*3 — 1, 13,13 -|- 1}, which are 


K,„ 


= i _ 27^pf2 _ 27^pf3 _ 27^13 , ^ 

3 3 3 3/i2 3/i2 3 /i 2 /i2 ’ 


K I 7*^1 <?? 7Pl2‘?2 7 ^Pi2 7^ 7P13<?3 7^13 

6 ^ 3/1 3/l2 6 /i2 ^ 3/1 3/l2 ’ 

K _ I 7‘?2 ^2 1^1 7^Pi 2 7^ 7P23<?3 7^pi3 

*1,*2±1.*3 6 3/1 3 /r 2 6/i 2 + 3/i 3 /i 2 ’ 

K I 7‘?3 <?! 7P13 7^Pi3 7^ 7P23<?2 7^13 

*1,Z2,Z3±1 g^ 6 3/1 3/l2 6 /i2 + 3/1 3/l2 ’ 


^2 T ^1 , ^l‘?2 7 

—1,13 7 ^ 


12/i 12 12/i2 


7Pl2- 


'?1 T ^2 2 P12 T P12 ± P13/323 


6h 


7 


6/l2 


^2±^1 ^l<r2 7^ , <?1±^2 2 /O 12 ± P 12 T P13P23 

A:*i±i.i2+1.*3 =7—TTTT— T -TTT “ TTiT? + 7Pi2—77- 7 


12/1 ^ 12 12/i2 

^3 T <?! , ^1<?3 7^ 

ttilil,12,23 — 1 “ 7 TOJ. ^ 10 


6/1 


6/l2 


12h 12 12/i2 


7P13 


“^1 T ^3 2 Pl3 T Pl3 ± P 12 P 23 


6h 


7 


6/l2 


<?3±^1 ^l<r3 7^ , <?1±^3 2 / 0 l 3 ± Pl3 T P 12 P 23 

^ii±i.i2.i3+i =7 —TTTT— T -ttt “ TATA + ^P^^—TT -7 


12/1 12 12/i2 

^3 T ‘?2 , ‘^2‘?3 7^ 

ttii,i2±l,j3 —1 7 TOJ. ^ 10 


6/1 


6/l2 


I2h 12 12/i2 


7P23 


^2 T <?3 2 P23 T P23 ± P 12 P 13 


6h 


7 


6/l2 


r> ‘^3 ± ^2 <?2C3 7^ , ^2 ± ^3 2 P 23 ± ^23 T P 12 P 13 

^ii,i2±i.i3+i =7 T -TTT “ 77777 + -7 


12/1 12 12/i2 


6h 


6/l2 


7> , P237i + P1372 + P1273 2 P23 T Pl2 T P13 2 P12P13 T P12P23 T P13P23 

— ± 7 7 03 1.2 7 


24/i 


24/i2 


12/l2 


^ P237i + P1372 + P1273 , 2 P23 T Pl2 ± Pl3 , 2 Pl2Pl3 ± P 12 P 23 T P13P23 

A: 7 ±i,i 2 +I,i 3-1 =T7-771-^7 - 777:7 -^7 


24/i 


24/i2 


12/l2 


jy P237i + P1372 + P1273 , 2P23±Pl2TPl3 , 2 P13P23 T P 12 P 23 ± P 12 P 13 

■?^ii±i.i2-i.*3+i = T 7-771-+ 7 -7777^-+ 7 


24 /i 


24/i2 


12/l2 


I P2371 + P1372 + P1273 2P23±P12±P13 2 P12P23 ± Pl2Pl3 ± P13P23 

A7±l,i2-|-1,43-1-1 — T7 03 1. ~ 7 03 1.2 ~7 


24/i ' 24/i2 

Similarly, we get the coefficients Mk,i,m of drUk,i,m (r), given by 


12/l2 


^ ^ P\3 '' '' P‘23 


^ ^ P\2 


—10 ^ 10 -,’ 


1 hc,2 
12 ^ 1^’ 


,7 _ 1 h<,i 

“12 ^ 127’ 

~Y2 ^ 1^’ ■^ 71.”1 = 2> 
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— — 0 — 1 , 771—1 —+ —1 — 0 - 

Additionally, we have g{x,T) = 0. We obtain a semi-discrete scheme of the form (fT^ . where 
Kx and M^ are time-independent. As in [5] we apply Crank-Nicolson type time discretisation and 
obtain the fully discrete scheme for the spatial interior. 

8.4 Treatment of the boundary conditions 

After deriving a high-order compact scheme for the spatial interior we now discuss the boundary 
conditions. 


8.4.1 Lo-wer boundaries 

The first boundary we discuss is = 0 for some i S / C {1,... n} at time t G [0,T[. Once the 
value of the asset is zero, it stays constant over time, see (HU. Hence, using = 0 for i G / in 
and applying the transformation (ED leads to 

9 n r^9 n n 

7 O U 2 ^ 

^ ^'^dx,dxj 

1 — 1 ^ = 1 1—1 

i<j 

with / = —Ut- Hence, at these boundaries we are able to obtain high-order compact schemes in 
the same manner as shown for the spatial interior with then n — |/| spatial dimensions, as the 
coefficients of the partial differential equations of these boundaries satisfy condition dill) . The 
case I = {!,..., n}, i.e. |/| = n, leads to the Dirichlet boundary condition ..., t) = 

• ■ • ! time T g]0 , Tmax], since in that case Ut = 0 . 


8.4.2 Upper boundaries 

Upper boundaries are boundaries with Si = Sf"^^ for some i G J C {1,..., n} at time t G [0,T[. 
For a sufficiently large for i € J, we set 

dV{S,,...,Sn,t) 

f)Q. ’ 

CfkJ'i ^._^max 

with Sk G ['S'™™, for k = {1,..., n} \ {i} for a European Power Put Basket option. Em¬ 

ploying this in dl6l) and using the transformations (flTll . yields 


2 ! 

2 ^ dx'j 
1—1 ^ 
i^J 


n 


7 ^ Y. 

7,1 = 1 


d'^u 

dxidxj 


i<j 


n 

aw, 

i^J 


( 21 ) 


with / = —Ut- Hence the upper boundaries show the same behaviour as the lower bound¬ 
aries for a European Power Put Basket. Analogously, we have the Dirichlet boundary condition 
ll(^max 5 ■ • ■ 5 ^Imax, 1”) ~ ll(3^max5 - - - j ^max? 0) for T G]0, Tjjiax] if T = {1, . . . , Tl}. 


8.5 Combination of npper and lower bonndaries 

A combination of upper and lower boundaries thus behaves in the same manner and the resulting 
partial differential equations with n — \I\ — | J| spatial dimensions satisfy condition dill) as well. 
Eor the corner points of O we have |J| -j- \J\ = n and thus again u = uq. 
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9 Numerical experiments for Black-Scholes Basket options 


In this section we discuss the numerical experiments for the Black-Scholes Basket Power Puts in 
spatial dimensions n = 2,3. The equation systems which have to be solved over time have been 
derived in Section |S1 According to [13] , we cannot expect fourth-order convergence if the initial 
condition is not sufficiently smooth. Hence, we have to smooth the initial condition for Power 
Puts with p = 1,2. In [T3] suitable smoothing operators are identified in Fourier space. Since the 
order of convergence of our high-order compact scheme is four, we use the smoothing operator $ 4 , 
given by its Fourier transform 


$4(0;) = 


sin(w/ 2 ) 

co/2 


4 r 


1 -I- - sin^(a;/2) 
o 


This leads to the smoothed initial condition 

3h 3h 


Uo (xi,X 2 ) = 


1 


$4 (0 $4 (0 Uo (xi -X,X 2 - y) dec dy. 


-3h -3h 


in the case n = 2 for any step size h > 0, where uq is the original initial condition and $ 4 (x) 
denotes the Fourier inverse of $4(0;), see [13]. If uq is smooth enough in the integrated region 
around (eci,..., Xn), we have uq (xi,..., x„) = uq (xi,..., x„). That means that it is possible to 
identify the points where smoothing is necessary. 




Figure 1: Example of grid points selected for the smoothing procedure in two space dimensions. 
We employ the smoothing operators of Kreiss et al. m to ensure high-order convergence of the 
approximations of the smoothed problem to the true solution of (fT31) . 

Figure [T] shows an example of a two-dimensional grid on the left side and on the right side 
a graph of the non-differentiable points of the initial condition given in ( 1201 ) together with the 
identified grid points, where smoothing is necessary. The points are chosen in such a way that 
we ensure that the non-differentiable points have no influence on uq (xi, X 2 ) for those points, 
which are not shown in Figure |T] on the right hand side. This approach reduces the necessary 
calculations significantly. As > 0, the smooth initial condition uq converges towards the original 
initial condition uq given in (1201) . The results in |13] guarantee high-order convergence of the 
approximation of the smoothed problem to the true solution of (IlSp . 

We use the relative Z^-error 11C/ref — [/||; 2 /||C/ref||/ 2 , as well as the Z°°-error ||t7ref — U\\io^ to 
examine the numerical convergence rate, where Uret denotes a reference solution on a fine grid and 
U is the approximation. When identifying the convergence order of the schemes, we determine it 
as the slope of the linear least square fit of the individual error points in the loglog-plots of error 
versus number of grid points per spatial direction. 

9.1 Numerical example with two underlying assets 

In this section we report the numerical results for a two-dimensional Black-Scholes Basket Power 
Put. We compare the high-order compact scheme (‘HOC’) with the standard scheme (‘2nd order’). 
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which is obtained by using the central difference operator directly in m for n = 2 with no further 
action and thus leads to a classical second-order scheme. We consider plain European Puts {p = 1) 
and use the smoothing procedure outlined above for the initial condition (1201) . The parameter 
values 


0-1 = 0.25, (72 = 0.35, 7 = 0.25, r = ln(1.05), wi = 0.35 = 1 - a;2, K = 10, 

and (5i = (52 = 0 are used, unless stated otherwise. The parabolic mesh ratio is fixed to Ar/h'^ = 
0.4, although we point out that neither the von Neumann stability analysis nor our numerical 
experiments revealed any practical restrictions on its choice. 



Figure 2: l°°- (left) and relative Z^-error (right) for two-dimensional Black-Scholes Basket Put and 
smoothed initial condition. 


Figure [2] shows convergence plots for the Z°°-error (left) and for the relative Z^-error (right) for a 
European Put, respectively. The initial condition is smoothed using the procedure outlined above. 
For both types of errors we observe that the numerical convergence rates agree very well with the 
theoretical orders of the schemes. The high-order compact scheme yields numerical convergence 
orders close to four and strongly outperforms the standard second-order scheme. The choice of 
the correlation parameter pi 2 = —0.8, pi2 = 0 and pi 2 = 0.8 has very little influence. 


9.2 Numerical example with three assets 

In this section we report on numerical experiments with three underlying assets. We choose the 
parameters 

,5i = 0.01, cr, =0.3, Wi = l/3, r = ln(1.05), 7 = 0.3, T = 0.25, K = 10. 

Due to the computational intensity of the three-dimensional problem the number of grid points 
per spatial dimension is smaller compared to the results in two dimensions reported above. To 
ensure that at the same time there is a sufficiently large number of grid points in time, we fix 
the parabolic mesh ratio to Ar/h? =0.1 (not for stability reasons). We perform two types of 
experiments: without any correlation between the assets (labeled by ‘nc’ in the plots), and with 
correlation (labeled by ‘c’ in the plots) using the parameter values pi,2 = —0.4, pi_3 = —0.1, 
P2,3 = —0.2. 

We compare the standard approximation to our high-order compact scheme for European 
Power Put options with p = 3,4. For the European Power Puts with p = 1, 2 one would smooth 
the initial condition, similar as above, to ensure high-order convergence. Figure [3] shows the 
convergence of the relative Z^-error for a European Power Put with p = 3 and p = 4. We use the 
original initial conditions, no smoothing is applied here. The numerical convergence rates of the 
high-order compact scheme are slightly reduced to about three and three and a half, respectively. 
Additional smoothing, which we omitted here due to limit the computational load, would result in 
even better results. Still, in the high-order compact scheme outperforms the standard second-order 
scheme significantly in all cases. 
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Figure 3: Relative ^^-error for three-dimensional Black-Scholes Basket Power Put, with p = 3 
(left) and p = 4 (right) 


9.3 Numerical example with space-dependent coefficients 

In this section we will apply numerical examples for where the continuous dividends are 

dependent on the underlying asset price. For both asset prices Si with i = 1,2 we consider 
the following example, where the continuous dividends are zero for small asset prices and then 
smoothly increase around an asset price S* > 0 towards a given parameter S* > 0, 

S- = S (S ) = ~ ~ 

Financially, the interpretation could be as follows: if the asset is a dividend-paying stock, low 
stock prices may mean that the company may not be in the financial position to pay dividends. A 
low value of Ci > 0 leads to slow transition from 0 to d*. We can apply the transformations given 
in and hence use the coefficients 


I J, 2 

— 2 ’ “ 7 Piji 


c.=7 y 


CTi r - 6i{Ke ') 


5 = 0 , 


( 22 ) 


for i = 1,2 to obtain the coefficients of the numerical scheme, see Section [5.11 The boundary 
conditions of Section 18.41 are employed and the parameter values of Section 19.11 as well as 


SI = 0.02, = 0.01, Cl = 0.35, C 2 = 0.5, S* = 0.9K/uj„ 


for i = 1,2 are used in the numerical experiments. Figure 0] shows numerical convergence plots 




Figure 4: l°°- (left) and relative Z^-error (right) for two-dimensional Black-Scholes Basket Put with 
space-dependent dividend and smoothed initial condition. 

for a European Put with space-dependent continuous dividend. Again, smoothing of the initial 
condition is employed. For the /°°-error as well as the Z^-error the high-order compact scheme 
has convergence rates close to four for pi _2 = 0, and pi _2 = 0.8. The convergence rate for the 
case pi ,2 = —0.8 is 3.22 in the Z°“-error, which is mainly due to the two approximations with 
eleven and 21 grid-points per spatial direction, and 3.57 in the Z^-error. The convergence orders 
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of the standard scheme are for pi ,2 = 0,0.8 are slightly above two for both types of errors. For 
pi ^2 = —0.8 the convergence orders are noticeable lower as well. In all cases of correlation the 
high-order compact scheme outperforms the standard second-order scheme significantly. 

10 Conclusion 

We presented a new high-order compact scheme for a class of parabolic partial differential equations 
with time and space dependent coefficients, including mixed second-order derivative terms in n 
spatial dimensions. The resulting schemes are fourth-order accurate in space and second-order 
accurate in time. In a thorough von Neumann stability analysis, where we focussed on the case of 
vanishing mixed derivative terms, we showed that a necessary stability condition holds for frozen 
coefficients without further conditions in two and three space dimensions. For non-vanishing mixed 
derivative terms we were able to give partial results. The results suggest unconditional stability of 
the scheme. As an application example we considered the pricing of European Power Put options 
in the multidimensional Black-Scholes model. The typical initial conditions of this problem lack 
sufficient regularity, therefore a suitable smoothing procedure was employed to ensure high-order 
convergence. In all numerical experiments performed a comparative standard second-order scheme 
is significantly outperformed. 

Although we derived the scheme in arbitrary space dimension, it was not our aim in this paper 
to attack the so-called curse of dimensionality. The issue of exponentially increasing number of 
unknowns with growing spatial dimension on full grids is of course alleviated to some degree by 
a high-order scheme. To obtain a similar accuracy as a second-order scheme which uses 0{N‘^) 
unknowns on a full grid, our high-order compact approach will ‘only’ require unknowns. 

To really attack very high-dimensional problems one would need to combine our approach with 
hierarchical approaches, e.g. using sparse grids (typically requiring 0(A^ln(A^)‘‘*“^) unknowns), 
which is beyond the scope of the present paper. 
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Appendix A Coefficients for semi-discrete scheme in three 
dimensions 


Considering an interior grid point S and time t G Or, the coefficients Kk,i,m 

of Uk,i,m (t) for A: e {zi — 1, fi, + 1}, I G {z2 — 1, *2, *2 + 1} and m G {13 — 1, is, is + 1} of the 
three-dimensional semi-discrete scheme in Section 5.2 are given by: 
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612 [a)26i3 [a)36i2 6i3[a)i6i2 [a)i623 ^ 612)623)2 

2462 T 48a6 ^ 48a6 ^ 48a2 6 ^ 48a6 96a6 

623)612)2 ^ C3612 ^ 613C2 ^ 613623 

96a6 48a6 48a6 24a62 ’ 


[c3)i 613 C3 [613)1 623[a)26i3 [613)11 [613)22 [613)33 

^ ^ “ m “ 126 ^ 24a2 6 ^ 48 ^ 48 ^ 48 


126 

^ L J-L ^ 

24a ^ 

48a ^ 48a 

24a 

^ 48a 662 

[a) 3 Cl 

_ 6i3[ci)i _ 

Ci[6i3)i 623 ) 613)23 1 

H 1 C 3 

C 1 C 3 H 2 [ 613)2 

24a 

^ 48a ^ 

48a 48a 

24a 

^ 24a 24a 

C 3[613 

[3 ^ 612 ) 613)12 

^ 613 ) 613)1 ^ 623 ) 613)2 

_6?3 

[a)i , C 3613 _ 6i3[a)3 

48a 

^ 48a 

24a6 24a6 

' 24a2/i 12ah ' 12ah 

612 [a); 

2613 ^ Cl 

613 ^ 613 ) 613)12 623 ) 01)2 

^ 623612 [a) 36 i 3 


48a 
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23 


fcisMiCi 

48a2 




Ei±l,j + l,m—l — 3“ 




bi2[o\i\biz\2 

48a2 


I ^23 [a]2[?)i3]3 
48a2 480^ 


&23 [a] 2 Cl 
48a2 

24 ■ 


, fcl3Ml[&13]3 
48a2 

&12H2[fel3]l 
48a2 ^ 


I &23[tt]3[^13]2 

48a2 

fclsMl _ Clbl 3 

V2ah 12ah 


fel2[ffl]2C3 fel3[fal3]3 

48a2 2Iah 

[C3]2 623 [&23]2 

^ ^ “ 12h 


6i2[C3]i ^ [a]3[623]3 &23[C2]2 

48a 24a ^ 48a ^ 

C3[&23]3 612 [ 623]12 &13[C2]l 

48a ^ 48a ^ 48a " 

612 H 1&23 ^ C 2 a ^ 623 C 3 
2Aa^h I2h 6h^ I2ah 

H2623 623[&23]23 C2[fe23]2 

24a2/i ^ 48a ^ 48a ^ 


fel3[a]3C3 b\2\bl3\2 

48a2 2Aah ’ 

C3 6 l 3 [a]lfe 23 [& 23 ]ll 

m ^ 24a2/i ^ 48 ■ 

&23[C2]2 &23[C3]3 H2C3 

48a ^ 48a 24a 


[623]22 
48 ^ 

[a]i[fe23]i 
24a 


[^ 23]33 
48 
C2C3 
^ 24a 


3 |_C 2 jl 6 i 3 [ 623 ]i 2 ^ [aj 2 l 023 j 2 ^ [aj 3 C 2 

48a ^ 48a 24a 24a 

fe23C3 & 23 M 3 &13[^23]l fc23[&23]2 &23 

12 a/i ^ I2ah 2Aah 2Aah ^ 12ah'^ 

b23]2 Cl [ 623 ]! ^ 612 &I 3 C 2&23 ^ 623 H 3 C 3 

I 8 a ^ 48a 12ah'^ ~ 12ah A8a^ 

6 l 2 [a]lC 3 6l3[a]l[&23]3 I &23[a]2[623]3 

48a2 48a2 48a2 

I &I 3 M 1 C 2 I fcl3[a]3[fc23]l I bl2[a]2[b23]l 

48a2 48a2 48a2 


[a] 2 [623)2 


623 [a] 3 [ 623)2 623 [a )2 , 

48a2 12 a 6 

623 M 2 C 2 , 612 [a)i [ 623)2 
48a2 48a2 


623)623)3 _ 612)623)1 6^3 [g)3 [623)3 [ 

24a6 24a6 2Aa‘^h 12h 

[613)2 613 [612)3 [623)1 612)613)1 

486 ^ 2462 T 486 ^ 486 ^ 96ah 

623612 6l3[a)3623 ^ 612^2623 623 

24a62 48a26 48a26 2462 ^ 


612 , 

2462 
623(612)2 
96a6 
[613)2 , 


a )2613 ^ 
48a6 
C 3612 
^ 48a6 


[0)3612 

48a6 
613C2 
^ 48a6 


613)0)1612 

48o2 6 
613623 
' 24o62’ 


[C2]3 

24 ’ 

613[612)1 623)613)3 

^ 96a6 ^ 96o6 

613)623)3 623C1 612613 

96a6 ^ 48o6 24a62 

1)1623 612(623)2 

48o6 ^ 96o6 


613)2 , 613 

486 2462 


623612 I ' 
24o62 
^12 I 
2462 z 
623(612)2 
96o6 

Mi ± h 

24 6 / 

[ 613)3 [ 

126 

[0)3 Cl ? 

24o 

C3[613)3 , 
48a 

612[g)2613 
24a2 6 
613M1C1 
48a2 


613 [612)3 [623)1 

2462 T 4g^ =F 
613)0)3623 I 612)0)2623 
48a26 48a26 

0)2613 [0)3612 613) 


48a6 4 

! C 3612 

^ 48a6 ^ 
613 C3 
662 + 126 ’ 
[0)3)613)3 , 
24a 

^ 6i3[ci)i 

' 48a 

I 612 ) 613)12 

48a 


0)3612 613)0)1612 


'12)613)1 613(612)1 623)613)3 

96a6 ^ 96a6 ^ 96a6 

623 613)623)3 623C1 612613 

2462 T 96a6 ^ 48a6 24a62 

'12 I [0)1623 612)623)2 


48a6 

613C2 

48a6 


48a2 6 
613623 
24a62’ 


[613)1 

.^623 [0)2613 

, [613)11 

, [613)22 

-H 

126 

^ 240^6 

48 

48 

48 

612)03)2 

, 613)03)3 _ 

[0)1)613)1 

I 02)613) 


48 a 

^ 48 o ^ 

24 a 

48 a 

662 


Cl[613)1 

48a 


623 ( 613)23 

48a 


[0)2)613)2 

24a 


613)613)1 

24a6 


623(613)2 

24a6 


613 ) 0)1 
24a2 6 


613 ) 0)3 

12 a 6 


613 613(613)12 623)01)2 623612 [0)3613 


126 12a62 48a 

612)0)1)613)2 623 [0)201 

48a2 ^ 48a2 


613 ) 0 ) 3 ) 613)1 

48a2 


612)0)203 

48a2 

623)0)203 

12 a 2 


613 ) 613)3 

24a6 

613 ) 0)103 

12 a 2 


623)0)2)613)3 
48a2 ^ 

_ 613)0)303 ^ 

48a2 

L [^3)3 [623)2 

6 66 


48a ' ; 

613)0)1)613)3 
' ^ 48a2 ^ 

612)0)2)613)1 

^ 48a2 

612)613)2 C1613 

24a6 12a6’ 


12a62 240^6 

623[0)3)613)2 
48a2 
613)0)1 
12 a 6 


6)03)11 

24 


6)03)33 

24 
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24 


± 

T 

T 

+ 

+ 

± 

± 

± 

KiJ±l,m+l = ± 
± 
± 


T 

T 

T 

+ 


h[c3]22 C 3 [a]ii [a ]22 [a ]33 ^ & 12 H 1623 ^ 612 H 2&13 

24 ^ m 12 12a2/i 12a2/i 

fe6l3Hl[C3]3 ft&12[a]l[C3]2 fe623H2[C3]3 ^^23 [«]3 [C3]2 ^IsHsHl 

24a2 ^ 24a2 240^ ^ 240^ 60 ^ 

fefci 3 M 3 [c 3 ]i hbi2[a]2[c3]i a &23 M 3 [a] 2 6 i 2 [a]i[a ]2 fei 3 [c 3 ]i 

24a2 ^ 24a2 ^ 60 ^ 60 ^ ^ 12 a 

Cl [g] 1 fejs ^12[a]l2 C 2 [a]2 &13 C3[a]3 bl3Ml2 I /i-C2[C3]2 

12a “ 6 a/i 2 ^ 12a 12a “ 6ah^ ~ 12a 12a 24a 

fc23H23 &23[C3]2 C2b23 , 623 [a] 2 &23[fe23]3 bl2[&23]l , fejsMs 

12 a 12 a ^ Qah 6ah ^ I2ah ^ 12ah I2a^h 

Ms^13 fclsMl fcl3[fal3]3 bl2[bl3]2 C 1613 ^ fc&13[c3]l2 hci [ 03)1 

I2a'^h Qah ^ 12ah ^ I2ah ^ Qah 24a 24a 

fe6l2[C3]l2 I hb23[C3]23 ^MlMll ftM2[C3l2 /i-MsMls 

24a 24a ^ 120 ^ l2a ^ 120 

feC3[C3]3 Ml Mi Mi 

24a 6 a 6 a 6 a ’ 

[C3]2 ^ &23 M3]2 C 3 6i3Mi^ 23 ^ [&23]ll ^ [623]22 [&23]33 

^ 12h m ^ 24a2/i 48 48 48 

6i2[C3]i MsMsJs &23[C2]2 ^ &23[C3]3 M 2 C 3 MlMsJl ^ C 2 C 3 

48a ^ Mo 48a 48a ^ 24a ^ Mo ^ 

CsMsJs ^ 6l2[fe23]l2 ^IsMll ^ 6 i3[623]i2 M 2 M 3 I 2 M 3 C 2 

48a 48a 48a 48a ^ 24a ^ 24a 

fel2Ml^23 I C2 g , fe23C3 & 23 M 3 1 &13[^23]l , fc23[fc23]2 02^23 

24a^h I2h Qh^ I2ah ^ I2ah 24ah 24ah 12ah 

M2fe23 &23 I fcssMsJsS , C2[fc23]2 , Cl [ 623 ] 1 &12bl3 & 23 M 2 

24a2/i I2ah^ 48a 48a 48a ^ 12ah'^ 12ah 

fc23M3M3]2 &23 M 3 C 3 &12Ml'^3 fclsMlMsjs ^23 [tt]2 [?)23]3 

isF ^ 48a2 ^ 48a2 ^ ^ 

& 23 M 2 C 2 6i2MiM3]2 &13Mi’^2 6i3M3[623]i &12M2[^23]i 

48a2 ^ isF ^ 48a2 ^ ^ 

^23[fc23]3 bl2[fc23]l _ ^23 M3 ^ [^ssjs ^ MJs 
24ah 24ah 24a^h 12h 24 


Note that in the above a, & 12 , 613 , & 23 , ci, C 2 , C 3 and g are evaluated at and 

r G Or. To streamline the notation we used [-Jfc and [•]kp to denote the first and second derivative 
of the coefficients with respect to Xfc, and with respect to Xk and Xp, respectively. 




